


 % Related parameters
aux  = (1-alpha)/((sigma^2)/2);
saux = s*lambda*aux;
 

 % Solve for unknown constant that first emerges in 
 % equ for density of surplus in natural wastage region. 
 % Constants in other regions linked to this unknown. 
 % Integrate density for array of possible values 
 % for this unknown. Solution is value that
 % ensure density integrates to 1.
 %
 % To come up with a grid of constants, use
 % distribution function from natural 
 % wastage region
%G_mM  = (A / saux) *(mL^aux2) *((mM/mL)^aux2 -1). 
%G_mM  < 1 & G_mM >0 -->
Ahi  = .99*saux*(mL^-saux)/((mM/mL)^saux -1);
Alo  = .01*saux*(mL^-saux)/((mM/mL)^saux -1);
obs_const= 500; 
constvec = linspace(Alo, Ahi, obs_const)';
 %
 % 
 % Elsby/Gottfries show that unknown constant shows up
 % in BC b/c it enters the density of surplus @ m=mL.
 % For any lambda and constant, U share of searchers
 % can be set to ensure BC is satisfied:
ushr  = ( (1/saux)*constvec*(mL^saux) ) ./ (1 + (1/saux)*constvec*(mL^saux));
 %
 %
 % Natural wastage (l,m)
obs_nw  = 300; 
[grid_nw, wghts_nw] = legendre(mL, mM, obs_nw);
mat_nw = grid_nw *ones(1,obs_const);
const = ones(obs_nw,1)*constvec'; 
g_nw = const.*(mat_nw.^(saux-1));   %obs_nw by obs_const
integral_nw = g_nw'*wghts_nw;       %obs_const by 1
G_mM = (constvec / saux) *(mL^saux) *((mM/mL)^saux -1);	%integral_nw;
 % J(m) in NW region:
J_nw = -(omega0/rho0) + ((1-omega1)/rho1)*grid_nw + J1*(grid_nw.^gamma1) + J2*(grid_nw.^gamma2);
 %     
 %
 %
 % Net Expansion (m,u)
obs_ne  = 200; 
[grid_ne, wghts_ne] = legendre(mM, mU, obs_ne);
 %
part1 = ((1-omega1)/(alpha*c))*(grid_ne - mM);
part2 = (r+(omega0/c))*(log(grid_ne) - log(mM));
delta_mM = s*lambda;
auxconst = ( ((1-omega1)/(alpha*c))*mM - (r+delta_mM + (omega0/c)) )*(mM^(-1/(1-alpha)));
part3 = auxconst*(1-alpha)*(grid_ne.^(1/(1-alpha)) - mM^(1/(1-alpha)));
 %
intdeltaoverm = (part1 - part2 - part3)*ones(1,obs_const);
G_ne = (ones(obs_ne,1)*G_mM').*exp(aux*intdeltaoverm) + ...
       (ones(obs_ne,1)*(ushr./(1-ushr))').*(exp(aux*intdeltaoverm) - 1); 
 %
gconst = G_mM + (ushr./(1-ushr));      
gridmat_ne = grid_ne*ones(1,obs_const);
deltamat_ne = ((1-omega1)/(alpha*c))*gridmat_ne - (r+(omega0/c)) ...
             - auxconst*(gridmat_ne.^(1/(1-alpha)));
g_ne = (ones(obs_ne,1)*gconst') .* aux.*(deltamat_ne./gridmat_ne) .* exp(aux*intdeltaoverm);     
 %
 %
 %
 % Sum up over two regions:
integral_of_g = G_ne(end,:)';
if min(integral_of_g) >1 || max(integral_of_g) < 1
    disp('ERROR: density cannot integrate to 1');
    pause;
end
 % Solution for unknown constant and related variables
constsol = interp1(integral_of_g, constvec, 1);
ushr     = ( (1/saux)*constsol*(mL^saux) ) ./ (1 + (1/saux)*constsol*(mL^saux));
urate    = s*ushr/(1 - ushr*(1-s));
L = En_firm/(1-urate);
 %
 %
 % Final density
g_nw = constsol*(grid_nw.^(saux-1));
g_ne = interp1(constvec, g_ne', constsol);	
g_ne = g_ne';
 % Final c.d.f.
G_mM = (constsol/ saux) *(mL^saux) *((mM/mL)^saux -1);	%integral_nw;
G_nw = cumsum( wghts_nw.*g_nw );                % (constsol / saux) *(mL^saux) *((grid_nw/mL).^saux -1);	% 
G_ne = G_mM + cumsum( wghts_ne.*g_ne);          % interp1(constvec, G_ne', constsol);                       % 
% G_ne = G_ne';
 % 
wghts_all= [wghts_nw;  wghts_ne];
grid_all = [grid_nw; grid_ne];
g_all = [g_nw; g_ne];
G_all = [G_nw; G_ne];
